TIP24 Analysis
Zach Peagler
2024-11-20
Introduction
This analysis has been performed to accompany a tomato growth trial in the Hydroponic Lab at the KSU Field Station (coordinates?). The trial was conducted with the purpose of determining the effects of differing soil microbial inoculation timings on tomato plant photosynthesis and crop yield and quality. Photosynthesis metrics were taken with a Li-COR Li-600 and two PhotosynQ MultispeQ v2.0s, and include stomatal conductance (gsw), measured in (mol m-2 s-1), which is a measure of how well the leaf is able to conduct gas through its stomates, and the efficiency of photosystem II (PhiPS2), where photosystem II is a piece of machinery in the chloroplast essential to photosynthesis and PhiPS2 is the quantum yield in light calculated from fluorescence.
Sample size (n=48) with 4 groups, 12 replicates per group. The tomato plants were grown in a dutch bucket hydroponic system from April to October 2024. Fluorometric measurements were taken twice weekly, and fruit were harvested upong fully ripening and taken back to the lab on the KSU campus for analysis.
Fruit analysis included weighing the tomato for its mass in grams and visually assessing it for blossom end rot (BER), a nutritional disorder caused by calcium deficiency that renders the fruit unmarketable, as well as visual checks for fungus and cracking. Fungus commonly appears on parts of the tomato with blossom end rot or cracking. From there, fruit were analyzed with a penetrometer for ripeness, (1-kg) where kg is mapped from 0:1, and then with Fisher Brix Refractometers for sugar concentration in percent sugar content of the tomato’s juice.
Discarded fruit were then composted.
1. Setup
1.1 Load packages
DATA WRANGLING - tidyverse
STATISTICS - lme4 - car - MASS - rstatix - pwr
GRAPHING - showtext - scico - ggpubr
OTHER - devtools - multiFitR - multiFitRgg - ztils
1.2 Graphical Setup
Load colors. Set the palette to one of the scientific color map palettes (scico) and create variants that contain the specific values for 2, 4, and 10 colors. This helps when working with factors, as it can be a pain to apply gradients otherwise. I also make the vector of colors =n+1 because the final value in scico palettes are often white, rendering them hard to see against a white plot background. (and we use white plot backgrounds because colored backgrounds are difficult to parse and expensive to print)
Load fonts. Google’s Open Sans for titles and Montserrat for body text.
Set shape selection from ggplot’s shape library.
a_palette <- "bilbao"
# super lazy way to make n+1 palettes
two_colors = scico(3, palette=a_palette)
four_colors = scico(5, palette=a_palette)
five_colors = scico(6, palette=a_palette)
true_two_col = scico(2, palette=a_palette)
ten_col = scico(10, palette=a_palette)
twelve_col = scico(13, palette=a_palette)
# add fonts
font_add_google("Open Sans", family = "open")
font_add_google("Montserrat", family = "mont")
# necessary to initialize imported fonts
showtext_auto()
# custom shapes
four_shapes = c(15,16,17,23)1.3 Experimental Constants
## germination date
germdate <- "2024-05-01"
## treatment order
treatment_order <- c("Control",
"Transplantation",
"Germination",
"Germ+Trans")
## plant factor order
plant_order <- c("1 1", "1 2", "1 3", "1 4", "1 5", "1 6",
"1 7", "1 8", "1 9", "1 10", "1 11", "1 12",
"2 1", "2 2", "2 3", "2 4", "2 5", "2 6",
"2 7", "2 8", "2 9", "2 10", "2 11", "2 12",
"3 1", "3 2", "3 3", "3 4", "3 5", "3 6",
"3 7", "3 8", "3 9", "3 10", "3 11", "3 12",
"4 1", "4 2", "4 3", "4 4", "4 5", "4 6",
"4 7", "4 8", "4 9", "4 10", "4 11", "4 12")2. Load Data
2.1 Photosynthesis data
2.1.1 Li-600 data
Load Li-600 data from file and convert strings to factors. Then change the treatments from numbers to the names of the groups. From here, filter the data for only those readings that have a leak percent of less than 10. QUESTION: What leak percent is acceptable to filter out? Rename the Date column to Date_ref to use it as a check that we formatted dates correctly. Format date to POSIXct in month-day-year format. Create a new plant factor that corresponds to each unique plant, then re-level the Treatment for consistency across analyses. Then we make a new column that is the log of gsw.
From here, outliers in the main target variables (stomatal conductance [gsw] and photosystem II efficiency [PhiPS2]) were identified.
This is one of those things I might not bother with, as it is almost certainly not necessary to remove outliers.
Then, new dataframes were created out of the base Li-600 dataset to look at the mean and sd of gsw and PhiPS2, grouped by 1. treatment and date and 2. plant factor.
Another dataframe was created looking at mass and fruit count sums by treatment and plant.
As a final note, this data has something wrong with it, as there are ~100 gsw values <0, which should be impossible. We’ve contacted Li-COR for clarification and hopefully there’s just a correction we can apply to the data. If not, then Li-COR has some explaining to do.
## Load and clean Li-600 data
Li_data_file <- "C:/Github/Portfolio/_data/tomato_inoculants/TIP24/TIP24_LI600.csv"
Li_data <- read.csv(Li_data_file, stringsAsFactors = T) %>%
mutate(Treatment = case_when(
Row==1~"Control",
Row==2~"Transplantation",
Row==3~"Germination",
Row==4~"Germ+Trans",
TRUE~NA)) %>%
filter(leak_pct<10, gsw>0) %>%
rename(Date_ref = Date) %>%
mutate(Date = parse_date_time(Date_ref, orders = "mdy"),
Time = parse_date_time(Time, orders = "T"),
DaysFromGermination = as.numeric(round(difftime(Date, germdate, units = c("days")), 0)),
plant_fac = as.factor(paste(Row, Pot)),
Treatment = factor(Treatment, levels = treatment_order),
logitPS2 = logit(PhiPS2, FALSE),
plant_fac = factor(plant_fac, levels = plant_order)
) %>%
group_by(DaysFromGermination) %>%
mutate(MinutesFromStart = round(difftime(Time, min(Time), units = "mins"), 4)) %>%
ungroup() %>%
mutate(
Time = format(Time, "%H:%M:%S")
)
# New dataframes of Li_data mean and sd for stomatal conductance (gsw) and Photosystem II efficiency (PhiPS2), grouped by treatment and plant
Li_data_stats_byplant <- Li_data %>%
group_by(Treatment, plant_fac) %>%
summarise_at(vars(gsw, logitPS2), list(mean=mean, sd=sd))2.1.2 MultispeQ data
mq_file <- "C:/Github/Portfolio/_data/tomato_inoculants/TIP24/TIP24_Multispeq.csv"
m_data <- read.csv(mq_file) %>%
mutate(Treatment = case_when(
Row=="A"~"Control",
Row=="B"~"Transplantation",
Row=="C"~"Germination",
Row=="D"~"Germ+Trans",
TRUE~NA),
Row_num = case_when(
Row=="A"~ 1,
Row=="B"~ 2,
Row=="C"~ 3,
Row=="D"~ 4,
TRUE~NA),
plant_fac = as.factor(paste(Row_num, Pot)),
Date = as.Date(time, "%m/%d/%Y"),
Time = parse_date_time(time, "%m/%d/%Y %H:%M %p"),
DaysFromGermination = as.numeric(round(difftime(Date, germdate, units = c("days")), 0)),
Treatment = factor(Treatment, levels = treatment_order),
plant_fac = factor(plant_fac, levels = plant_order),
logitPS2 = logit(Phi2, FALSE),
logitFvPFmP = logit(FvP_over_FmP, FALSE)
) %>%
group_by(DaysFromGermination) %>%
mutate(MinutesFromStart = round(difftime(Time, min(Time), units = "mins"), 4)) %>%
ungroup() %>%
mutate(
Time = format(Time, "%H:%M:%S")
)
m_data_stats_byplant <- m_data %>%
group_by(Treatment, plant_fac) %>%
summarise_at(vars(logitFvPFmP, logitPS2), list(mean=mean, sd=sd))2.2 Fruit Data
Load fruit data from set file name and location, converting strings to factors and renaming Treatment’s group names. Then setup “fruit” as a way to sum the total fruit count, create references for our dates then parse them to POSIXct in month-day-year format, then grab just the days and calculate the difference in days from harvest to analysis. (NOTE: this is susceptible to month rollover errors, so be sure not to make any of those mistakes) Make a plant factor as a unique ID for each plant, relevel the treatment for consistency across analyses, create mass bins (n=10) to look at BER occurrence across mass groups, and convert BER, fungus, and cracking to factors. Then, we map the penetrometer data from 0:1 and subtract it from 1.0, effectively giving us the ripeness. (this is because higher penetrometer values correspond to lower ripeness)
Then create a copy of the data that doesn’t have BER and another copy of just those with BER.
Make a pivot table of the means and sd of the tomato mass (g) and sugar (%), grouped by treatment.
Make another pivot table to get the sum of fruit count, mass, and BER, fungus, and cracking occurrence by treatment group. Then calculate the probability of BER, fungus, and cracking for each group. Do that again grouping the fruit into mass bins.
## Load and clean TIP24 fruit data
Fl_data_file <- "C:/Github/Portfolio/_data/tomato_inoculants/TIP24/TIP24_Fruit.csv"
Fl_data <- read.csv(Fl_data_file, stringsAsFactors = T) %>%
filter(mass >0) %>%
mutate(Treatment = case_when(
row==1~"Control",
row==2~"Transplantation",
row==3~"Germination",
row==4~"Germ+Trans",
TRUE~NA),
fruit = 1,
date_analysis_ref = date_analysis,
date_harvest_ref = date_harvest,
date_analysis = parse_date_time(date_analysis_ref, orders = "mdy"),
date_harvest = parse_date_time(date_harvest_ref, orders="mdy"),
DaysFromHarvestToAnalysis = as.numeric(round(difftime(date_analysis,
date_harvest, units = c("days")), 0)),
DaysFromGermination = as.numeric(round(difftime(date_analysis,
germdate, units = c("days")), 0)),
plant_fac = as.factor(paste(row, plant)),
Treatment = factor(Treatment, levels = treatment_order),
plant_fac = factor(plant_fac, levels = plant_order),
mass_bin = cut(mass, breaks=10),
BER_fac = as.factor(BER),
fungus_fac = as.factor(fungus),
cracking_fac = as.factor(cracking),
ripeness = abs(1 - round(penetrometer/max(na.omit(penetrometer)), 2))
)
## Filter to only rows where there is no BER and sugar is greater than 0
Fl_data_no_BER <- Fl_data %>%
filter(BER==0 & sugar_avg > 0)
## New dataframes, filtering to only those where BER is equal to 1
Fl_data_BER <- Fl_data %>%
filter(BER==1)
## Make summary dataframe with mean and sd of mass and sugar by plant.
Fl_means_byplant <- Fl_data_no_BER %>%
group_by(Treatment, plant_fac) %>%
summarise_at(vars(mass, sugar_avg), list(mean=mean, sd=sd))
## Sum the fruit, BER, fungus, cracking, and mass by treatment group, then get probs of each factor.
Fl_data_summary <- Fl_data %>%
group_by(Treatment) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking, mass),
list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4),
pfungus = round(fungus_sum/fruit_sum, 4),
pcracking = round(cracking_sum/fruit_sum, 4)
)
## Sum the fruit, BER, fungus, and cracking by mass bin, then get probs of each factor.
Fl_data_mb <- Fl_data %>%
group_by(Treatment, mass_bin) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking), list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4),
pfungus = round(fungus_sum/fruit_sum, 4),
pcracking = round(cracking_sum/fruit_sum, 4)
)
## Sum the fruit, BER, fungus, and cracking by plant, then get probs of each factor.
Fl_sum_byplant <- Fl_data %>%
group_by(Treatment, plant) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking), list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4),
pfungus = round(fungus_sum/fruit_sum, 4),
pcracking = round(cracking_sum/fruit_sum, 4)
)
Fl_sum_means <- Fl_sum_byplant %>%
group_by(Treatment) %>%
summarise_at(vars(fruit_sum, mass_sum, pBER), list(mean=mean, sd=sd))
Fl_means <- Fl_data %>%
group_by(Treatment) %>%
summarise_at(vars(mass), list(mean=mean, sd=sd))3. Analyses
3.1 Photosynthesis Analysis
3.1.1 Stomatal Conductance (GSW) exploratory
# PDF, CDF, and KS test
multiPDF_Z(Li_data$gsw, 100, "all", a_palette, "gsw")
multiCDF_Z(Li_data$gsw, 100, "all", a_palette, "gsw")
multiKS_cont(Li_data$gsw, "all")## Distribution Distance PValue
## 1 Normal 0.165 0.000
## 2 Lognormal 0.089 0.000
## 3 Gamma 0.058 0.039
## 4 Exponential 0.139 0.000# check for heteroscedasticity
## in gsw as a function of treatment
leveneTest(Li_data$gsw~Li_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.8428 0.03717 *
## 578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## in gsw as a function of time from start
leveneTest(Li_data$gsw~Li_data$MinutesFromStart)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 342 0.9056 0.7996
## 239## gsw by Treatment boxplot
ggplot(data = Li_data, aes(x= Treatment, y = gsw, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
geom_jitter( width=.15, height=0)+
# stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
# stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Stomatal Conductance Across Inoculation Treatments in Tomato", 40)
) +
ylab("Stomatal Conductance (mol m-2 s-1)")+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
# Mean GSW per plant by Treatment boxplot
ggplot(data = Li_data_stats_byplant, aes(x= Treatment, y = gsw_mean, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
geom_jitter( width=.15, height=0)+
# stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
# stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Mean Stomatal Conductance Per Plant
Across Inoculation Treatments in Tomato", 50)
) +
ylab("Stomatal Conductance (mol m-2 s-1)")+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
# GSW by date across relative humidity
ggplot(data = Li_data, aes(x=DaysFromGermination, y = gsw,
color = rh_s)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_scico(begin=0.9, end=0, palette=a_palette)+
# scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
theme_bw()+
ylab("Stomatal Conductance (mol m-2 s-1)")+
xlab("Days From Germination")+
labs(
title=str_wrap("Stomatal Conductance By Days From Germination Across Relative Humidity in Tomato", 50)
)+
guides(color=guide_colorbar(title=str_wrap("Relative Humidity", 10)))+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="right",
legend.title.position = "top",
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=two_colors[2], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
# GSW by time across relative humidity
ggplot(data = Li_data, aes(x=Time, y = gsw,
color = Date_ref,
shape = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_shape_manual(values=four_shapes)+
scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
theme_bw()+
ylab("Stomatal Conductance (mol m-2 s-1)")+
xlab("Time")+
labs(
title=str_wrap("Stomatal Conductance By Time Across Date in Tomato", 50)
)+
guides(color=guide_legend(title="Date"))+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="right",
legend.title.position = "top",
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=two_colors[2], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
# gsw by relative humidity across treatment groups
ggplot(data = Li_data, aes(x=rh_s, y = gsw,
shape=Treatment, fill=Treatment, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
#ylim(0,1)+
theme_bw()+
ylab("Stomatal Conductance (mol m-2 s-1)")+
xlab("Relative Humidity")+
labs(
title=str_wrap("Stomatal Conductance By Relative
Humidity Across Inoculation Treatment in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="inside",
legend.title.position = "top",
legend.position.inside=c(0.15,0.75),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
# gsw by Minutes from Start
ggplot(data = Li_data, aes(x=MinutesFromStart, y = gsw, color = rh_s)) +
geom_point()+
scale_color_scico(begin=0.9, end=0, palette=a_palette)+
facet_wrap(~Treatment)+
theme_bw()+
ylab("Stomatal Conductance (mol m-2 s-1)")+
xlab("Minutes From Start")+
labs(
title=str_wrap("Stomatal Conductance by Minutes From Start Across Inoculation Treatment in Tomato", 50)
)+
guides(color=guide_colorbar(title=str_wrap("Relative Humidity", 8)))+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.
3.1.3 Stomatal Conductance (GSW) Analysis
Do a generalized linear mixed model with stomatal conductance (gsw) as response, Treatment as predictor, and days from germination as a fixed effect that should incorporate fluctuations in environmental variables, such as relative humidity, with a gamma family, log link function because gsw best fits a gamma distribution.
#GLMM with Treatment as explanatory and days from germination as fixed effect
gsw_glmm_dfg_summary <- summary(gsw_glm_dfg <- (glm(gsw ~ Treatment + DaysFromGermination,
data=Li_data, family=Gamma(link="log"))))
print(gsw_glmm_dfg_summary)##
## Call:
## glm(formula = gsw ~ Treatment + DaysFromGermination, family = Gamma(link = "log"),
## data = Li_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.545480 0.120872 -4.513 7.75e-06 ***
## TreatmentTransplantation -0.215724 0.099400 -2.170 0.03039 *
## TreatmentGermination -0.285432 0.100647 -2.836 0.00473 **
## TreatmentGerm+Trans -0.541036 0.102431 -5.282 1.81e-07 ***
## DaysFromGermination 0.003901 0.002176 1.792 0.07361 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.7609672)
##
## Null deviance: 418.91 on 581 degrees of freedom
## Residual deviance: 395.12 on 577 degrees of freedom
## AIC: 396.87
##
## Number of Fisher Scoring iterations: 6#GLMM with Treatment as explanatory and relative humidity as fixed effect
gsw_glmm_rhs_summary <- summary(gsw_glm_rhs <- (glm(gsw ~ Treatment + rh_s,
data=Li_data, family=Gamma(link="log"))))
print(gsw_glmm_rhs_summary)##
## Call:
## glm(formula = gsw ~ Treatment + rh_s, family = Gamma(link = "log"),
## data = Li_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.903211 0.199662 -29.566 < 2e-16 ***
## TreatmentTransplantation -0.248352 0.057316 -4.333 1.73e-05 ***
## TreatmentGermination -0.269800 0.058169 -4.638 4.35e-06 ***
## TreatmentGerm+Trans -0.269977 0.059667 -4.525 7.35e-06 ***
## rh_s 0.071222 0.002588 27.519 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2532233)
##
## Null deviance: 418.91 on 581 degrees of freedom
## Residual deviance: 182.51 on 577 degrees of freedom
## AIC: -87.276
##
## Number of Fisher Scoring iterations: 9# pseudo R2 of glmms
gsw_glmm_dfg_pseudoR2 <- (gsw_glmm_dfg_summary$null.deviance -
gsw_glmm_dfg_summary$deviance)/gsw_glmm_dfg_summary$null.deviance
print(gsw_glmm_dfg_pseudoR2)## [1] 0.05678831gsw_glmm_rhs_pseudoR2 <- (gsw_glmm_rhs_summary$null.deviance -
gsw_glmm_rhs_summary$deviance)/gsw_glmm_rhs_summary$null.deviance
print(gsw_glmm_rhs_pseudoR2)## [1] 0.5643147# PREDICTIONS
## im going to do this in maybe the worst possible way, sorry
# GSW_RHS
## initialize prediction plot
p <- ggplot()+
geom_point(data=Li_data, aes(x=rh_s, y=gsw, color=Treatment, shape=Treatment),
alpha = 0.8)
## base GLM
gsw_rhs <- (glm(gsw ~ rh_s,
data=Li_data,
family=Gamma(link="log")))
## sequence
rh_s <- seq(min(Li_data$rh_s), max(Li_data$rh_s), length=200)
# loop over treatments
for (x in treatment_order) {
## treatment specific GLM
gsw_rhs <- (glm(gsw ~ rh_s,
data=Li_data %>% filter(Treatment == x), #treatmnent == x
family=Gamma(link="log")))
## calculate predictions
prx <- as.data.frame(x=rh_s)
pred <- predict(gsw_rhs, newdata=prx, se.fit=TRUE, type="link")
prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
prx$Treatment <- x
prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
## add predicted points
p <- p +
geom_line(data=prx, aes(x=rh_s, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_line(data=prx, aes(x=rh_s, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
geom_line(data=prx, aes(x=rh_s, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_ribbon(data=prx, aes(x=rh_s, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
facet_wrap(~Treatment)+
labs(title = "Real vs Predicted for gsw_rhs glm")+
ylab("Stomatal Conductance (mol+1 m-2 s-1)")+
xlab("Relative Humidity (%)")+
theme_minimal()+
theme(
text = element_text(size=24, family="mont"),
legend.position="none",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
print(p)
# GSW_DFG
## initialize prediction plot
p <- ggplot()+
geom_point(data=Li_data, aes(x=DaysFromGermination, y=gsw, color=Treatment, shape=Treatment),
alpha = 0.8)
## base GLM
gsw_dfg <- (glm(gsw ~ DaysFromGermination,
data=Li_data,
family=Gamma(link="log")))
## sequence
DaysFromGermination <- seq(min(Li_data$DaysFromGermination), max(Li_data$DaysFromGermination), length=200)
# loop over treatments
for (x in treatment_order) {
## treatment specific GLM
glm_x <- (glm(gsw ~ DaysFromGermination,
data=Li_data %>% filter(Treatment == x), #treatmnent == x
family=Gamma(link="log")))
## calculate predictions
prx <- as.data.frame(x=DaysFromGermination)
pred <- predict(glm_x, newdata=prx, se.fit=TRUE, type="link")
prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
prx$Treatment <- x
prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
## add predicted points
p <- p +
geom_line(data=prx, aes(x=DaysFromGermination, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_line(data=prx, aes(x=DaysFromGermination, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
geom_line(data=prx, aes(x=DaysFromGermination, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_ribbon(data=prx, aes(x=DaysFromGermination, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
facet_wrap(~Treatment)+
labs(title = "Real vs Predicted for gsw_dfg glm")+
ylab("Stomatal Conductance (mol+1 m-2 s-1)")+
xlab("Days From Germination")+
theme_minimal()+
theme(
text = element_text(size=24, family="mont"),
legend.position="none",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
print(p)
3.1.3 PhiPS2 exploratory
# Li-600
## PDF and CDF plots
multiPDF_Z(subset(Li_data$logitPS2, Li_data$logitPS2>0), 100, "all", a_palette, "Logit PhiPS2")
multiCDF_Z(subset(Li_data$logitPS2, Li_data$logitPS2>0), 100, "all", a_palette, "Logit PhiPS2")
## multi KS tests
multiKS_cont(subset(Li_data$logitPS2, Li_data$logitPS2>0), "all")## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "pexp", rate = x_exp$estimate): ties should not
## be present for the Kolmogorov-Smirnov test## Distribution Distance PValue
## 1 Normal 0.310 0
## 2 Lognormal 0.157 0
## 3 Gamma 0.212 0
## 4 Exponential 0.449 0## test for heteroscedasticity
leveneTest(Li_data$logitPS2~Li_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.036 0.1077
## 578leveneTest(Li_data$logitPS2~Li_data$MinutesFromStart)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 342 0.7906 0.9766
## 239# Li PhiPS2 by Days From Germination, colored by treatment
ggplot(data = Li_data, aes(x=DaysFromGermination, y = PhiPS2, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
ylim(.6,.8)+
theme_bw()+
ylab("PhiPS2")+
xlab("Days From Germination")+
labs(
title=str_wrap("Li-600 Photosystem II Efficiency by Date Across Inoculation Treatment in Tomato", 50)
)+
guides(color=guide_legend(title=str_wrap("Treatment", 8)))+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# PhiPS2 by Minutes from Start colored by treatment
ggplot(data = Li_data, aes(x=MinutesFromStart, y = PhiPS2, color = Treatment)) +
geom_point()+
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
theme_bw()+
ylab("PhiPS2")+
xlab("Minutes From Start")+
labs(
title=str_wrap("Li-600 Photosystem II Efficiency by Minutes From Start Across Inoculation Treatment in Tomato", 50)
)+
guides(color=guide_legend(title=str_wrap("Treatment", 8)))+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.
## logitPS2 by Treatment boxplot
ggplot(data = Li_data, aes(x= Treatment, y = logitPS2, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
geom_jitter( width=.05, height=0)+
# stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
# stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Li-600 LogitPS2 Across Inoculation Treatments in Tomato", 40)
) +
ylab("Logit PhiPS2")+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
# Mean logitPS2 per plant by Treatment boxplot
ggplot(data = Li_data_stats_byplant, aes(x= Treatment, y = logitPS2_mean, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
geom_jitter( width=.15, height=0)+
# stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
# stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Li-600 Mean Logit PhiPS2 Per Plant
Across Inoculation Treatments in Tomato", 60)
) +
ylab("Logit PhiPS2")+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
# MultispeQ
## PDF, CDF, and KS test
multiPDF_Z(subset(m_data$logitPS2, m_data$logitPS2>0),
100, c('normal','lognormal','gamma'), a_palette, "Logit PhiPS2")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
multiCDF_Z(subset(m_data$logitPS2, m_data$logitPS2>0),
100, c('normal','lognormal','gamma'), a_palette, "Logit PhiPS2")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
multiKS_cont(subset(m_data$logitPS2, m_data$logitPS2>0),
c('normal','lognormal','gamma'))## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Distribution Distance PValue
## 1 Normal 0.048 0.070
## 2 Lognormal 0.052 0.043
## 3 Gamma 0.034 0.393# check for heteroscedasticity
leveneTest(m_data$logitPS2~m_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 3.6607 0.01224 *
## 718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggplot(data = m_data, aes(x=DaysFromGermination, y = logitPS2,
shape=Treatment, fill=Treatment, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
theme_bw()+
ylab("Logit PhiPS2")+
xlab("Days From Germination")+
labs(
title=str_wrap("MultispeQ Logit Photosystem II Efficiency By Date Across
Inoculation Treatment in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
ggplot(data = m_data, aes(x=MinutesFromStart, y = logitPS2,
shape=Treatment, fill=Treatment, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
theme_bw()+
ylab("Logit PhiPS2")+
xlab("Minutes From Start")+
labs(
title=str_wrap("MultispeQ Logit Photosystem II Efficiency By Minutes From Start Across
Inoculation Treatment in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.
## logitPS2 by Treatment boxplot
ggplot(data = m_data, aes(x= Treatment, y = logitPS2, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
geom_jitter( width=.15, height=0)+
# stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
# stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("MultispeQ LogitPS2 Across Inoculation Treatments in Tomato", 40)
) +
ylab("Logit PhiPS2")+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
# Mean logitPS2 per plant by Treatment boxplot
ggplot(data = m_data_stats_byplant, aes(x= Treatment, y = logitPS2_mean, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.3, outliers = FALSE)+
geom_jitter( width=.15, height=0)+
# stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
# stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("MultispeQ Mean Logit PhiPS2 Per Plant
Across Inoculation Treatments in Tomato", 50)
) +
ylab("Mean Logit PhiPS2")+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## PhiPS2 by date, grouped by treatment
ggplot(data = Li_data, aes(x=as.factor(Date), y = PhiPS2, color = Treatment,
fill=Treatment, shape=Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
ylim(.6,.8)+
facet_wrap(~Treatment)+
theme_bw()+
ylab("PhiPS2")+
xlab("Date")+
labs(
title=str_wrap("PhiPS2 By Date Across Microbial Treatments in Tomato", 60)
)+
theme(
text = element_text(size=24, family="mont"),
legend.position="none",
axis.title = element_text(size=24, family = "mont", face= "bold"),
axis.text.x = element_text(size=12, family = "mont", angle=10),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# PhiPS2 by leaf across temperature
ggplot(data = Li_data, aes(x=Date, y = PhiPS2, color = Tleaf)) +
geom_jitter(width=1, size=2, alpha=1, aes(x=as.factor(Date)))+
scale_color_scico(begin=0.9, end=0, palette=a_palette)+
scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
ylim(.6,.8)+
theme_bw()+
guides(color=guide_legend(title="Leaf Temp (C)"))+
ylab("PhiPS2")+
xlab("Date")+
labs(
title=str_wrap("PhiPS2 By Date Across Temperature in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="inside",
legend.title.position = "top",
legend.position.inside=c(0.75,0.2),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=two_colors[2],
fill=NA, linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
3.1.4 PhiPS2 Analysis
Because we logit transformed PhiPS2, we should be able to treat it as normal.
# linear mixed model with PhiPS2 as response, Treatment as predictor (explanatory), days from germination as fixed effect
mod_Phi2_dfg_LI600_s <- summary(mod_Phi2_dfg_LI600 <- (lm(
logitPS2 ~ Treatment + DaysFromGermination,
data = Li_data)))
print(mod_Phi2_dfg_LI600_s)##
## Call:
## lm(formula = logitPS2 ~ Treatment + DaysFromGermination, data = Li_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9216 -0.1892 -0.0868 0.0226 4.4045
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.477569 0.079358 6.018 3.14e-09 ***
## TreatmentTransplantation 0.158414 0.065261 2.427 0.01551 *
## TreatmentGermination 0.097558 0.066080 1.476 0.14039
## TreatmentGerm+Trans 0.203130 0.067251 3.020 0.00264 **
## DaysFromGermination 0.008719 0.001429 6.102 1.93e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5727 on 577 degrees of freedom
## Multiple R-squared: 0.07766, Adjusted R-squared: 0.07126
## F-statistic: 12.15 on 4 and 577 DF, p-value: 1.741e-09## R2 0.0713
mod_Phi2_qamb_LI600_s <- summary(mod_Phi2_dfg_LI600 <- (lm(
logitPS2 ~ Treatment + Qamb,
data = Li_data)))
print(mod_Phi2_qamb_LI600_s)##
## Call:
## lm(formula = logitPS2 ~ Treatment + Qamb, data = Li_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8591 -0.1740 -0.0763 0.0093 4.6414
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0644152 0.0889739 11.963 < 2e-16 ***
## TreatmentTransplantation 0.2112193 0.0687409 3.073 0.00222 **
## TreatmentGermination 0.0784828 0.0688119 1.141 0.25453
## TreatmentGerm+Trans 0.2283420 0.0695834 3.282 0.00109 **
## Qamb -0.0006152 0.0002478 -2.483 0.01333 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5878 on 577 degrees of freedom
## Multiple R-squared: 0.02852, Adjusted R-squared: 0.02179
## F-statistic: 4.235 on 4 and 577 DF, p-value: 0.002186## R2 0.0218
mod_Phi2_dfg_Mult_s <- summary(mod_Phi2_dfg_Mult <- (lm(
logitPS2 ~ Treatment + DaysFromGermination,
data = m_data)))
print(mod_Phi2_dfg_Mult_s)##
## Call:
## lm(formula = logitPS2 ~ Treatment + DaysFromGermination, data = m_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60212 -0.10424 -0.00623 0.09618 0.46918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4217061 0.0228180 18.481 < 2e-16 ***
## TreatmentTransplantation 0.0968281 0.0155648 6.221 8.4e-10 ***
## TreatmentGermination -0.0394771 0.0154092 -2.562 0.0106 *
## TreatmentGerm+Trans 0.0394109 0.0156349 2.521 0.0119 *
## DaysFromGermination 0.0004075 0.0002951 1.381 0.1677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1488 on 717 degrees of freedom
## Multiple R-squared: 0.1048, Adjusted R-squared: 0.09982
## F-statistic: 20.99 on 4 and 717 DF, p-value: 2.223e-16## R2 0.0998
mod_Phi2_qamb_Mult_s <- summary(mod_Phi2_qamb_Mult <- (lm(
logitPS2 ~ Treatment + Light.Intensity..PAR.,
data = m_data)))
print(mod_Phi2_qamb_Mult_s)##
## Call:
## lm(formula = logitPS2 ~ Treatment + Light.Intensity..PAR., data = m_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40277 -0.06556 -0.00085 0.07001 0.32258
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.129e-01 1.344e-02 53.039 < 2e-16 ***
## TreatmentTransplantation 3.563e-02 1.179e-02 3.021 0.00261 **
## TreatmentGermination -8.851e-02 1.159e-02 -7.639 7e-14 ***
## TreatmentGerm+Trans -2.665e-02 1.189e-02 -2.241 0.02530 *
## Light.Intensity..PAR. -9.182e-04 3.767e-05 -24.377 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1102 on 717 degrees of freedom
## Multiple R-squared: 0.5092, Adjusted R-squared: 0.5065
## F-statistic: 186 on 4 and 717 DF, p-value: < 2.2e-16## R2 0.5065
# run and plot predictions
# Li-600 logitPS2 ~ Qamb
## initialize prediction plot
p <- ggplot()+
geom_point(data=Li_data, aes(x=Qamb, y=logitPS2, color=Treatment, shape=Treatment),
alpha = 0.8)
## base LM
lm_x <- (lm(logitPS2 ~ Qamb, data=Li_data))
## sequence
Qamb <- seq(min(Li_data$Qamb), max(Li_data$Qamb), length=200)
# loop over treatments
for (x in treatment_order) {
## treatment specific GLM
lm_x <- (lm(logitPS2 ~ Qamb,
data=Li_data %>% filter(Treatment == x)))
## calculate predictions
prx <- as.data.frame(x=Qamb)
pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
prx$Treatment <- x
prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
## add predicted points
p <- p +
geom_line(data=prx, aes(x=Qamb, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_line(data=prx, aes(x=Qamb, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
geom_line(data=prx, aes(x=Qamb, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_ribbon(data=prx, aes(x=Qamb, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for Li-600 logitPS2 ~ Treatment + Qamb")+
ylab("logit PhiPS2")+
xlab("Ambient Light (lumens)")+
theme_minimal()+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
print(p)
# Li-600 logitPS2 ~ DaysFromGermination
## initialize prediction plot
p <- ggplot()+
geom_point(data=Li_data, aes(x=DaysFromGermination, y=logitPS2, color=Treatment, shape=Treatment),
alpha = 0.8)
## base LM
lm_x <- (lm(logitPS2 ~ DaysFromGermination,
data=Li_data))
## sequence
DaysFromGermination <- seq(min(Li_data$DaysFromGermination), max(Li_data$DaysFromGermination), length=200)
# loop over treatments
for (x in treatment_order) {
## treatment specific GLM
lm_x <- (lm(logitPS2 ~ DaysFromGermination,
data=Li_data %>% filter(Treatment == x)))
## calculate predictions
prx <- as.data.frame(x=DaysFromGermination)
pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
prx$Treatment <- x
prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
## add predicted points
p <- p +
geom_line(data=prx, aes(x=DaysFromGermination, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_line(data=prx, aes(x=DaysFromGermination, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
geom_line(data=prx, aes(x=DaysFromGermination, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_ribbon(data=prx, aes(x=DaysFromGermination, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for Li-600 logitPS2 ~ Treatment + DaysFromGermination")+
ylab("logit PhiPS2")+
xlab("Days From Germination")+
theme_minimal()+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
print(p)
# Multispeq logitPS2 ~ DaysFromGermination
## initialize prediction plot
p <- ggplot()+
geom_point(data=m_data, aes(x=DaysFromGermination, y=logitPS2, color=Treatment, shape=Treatment),
alpha = 0.8)
## base GLM
lm_x <- (lm(logitPS2 ~ DaysFromGermination,
data=m_data))
## sequence
DaysFromGermination <- seq(min(m_data$DaysFromGermination), max(m_data$DaysFromGermination), length=200)
# loop over treatments
for (x in treatment_order) {
## treatment specific GLM
lm_x <- (lm(logitPS2 ~ DaysFromGermination,
data=m_data %>% filter(Treatment == x))) #treatmnent == x
## calculate predictions
prx <- as.data.frame(x=DaysFromGermination)
pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
prx$Treatment <- x
prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
## add predicted points
p <- p +
geom_line(data=prx, aes(x=DaysFromGermination, y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_line(data=prx, aes(x=DaysFromGermination, y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
geom_line(data=prx, aes(x=DaysFromGermination, y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_ribbon(data=prx, aes(x=DaysFromGermination, ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for MultispeQ logitPS2 ~ Treatment + DaysFromGermination")+
ylab("logit PhiPS2")+
xlab("Days From Germination")+
theme_minimal()+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
print(p)
# Multispeq logitPS2 ~ Light.Intensity..PAR.
## initialize prediction plot
p <- ggplot()+
geom_point(data=m_data, aes(x=Light.Intensity..PAR., y=logitPS2, color=Treatment, shape=Treatment),
alpha = 0.8)
## base GLM
lm_x <- (lm(logitPS2 ~ Light.Intensity..PAR.,
data=m_data))
## sequence
Light.Intensity..PAR. <- seq(min(m_data$Light.Intensity..PAR.), max(m_data$Light.Intensity..PAR.), length=200)
# loop over treatments
for (x in treatment_order) {
## treatment specific GLM
lm_x <- (lm(logitPS2 ~ Light.Intensity..PAR.,
data=m_data %>% filter(Treatment == x))) #treatmnent == x
## calculate predictions
prx <- as.data.frame(x=Light.Intensity..PAR.)
pred <- predict(lm_x, newdata=prx, se.fit=TRUE, type="terms")
prx$lo <- exp(qnorm(0.025, pred$fit, pred$se.fit))
prx$mn <- exp(qnorm(0.5, pred$fit, pred$se.fit))
prx$up <- exp(qnorm(0.975, pred$fit, pred$se.fit))
prx$Treatment <- x
prx$Treatment <- factor(prx$Treatment, levels = treatment_order)
## add predicted points
p <- p +
geom_line(data=prx, aes(x=Light.Intensity..PAR., y=lo, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_line(data=prx, aes(x=Light.Intensity..PAR., y=mn, color = Treatment), linewidth=2, show.legend=FALSE)+
geom_line(data=prx, aes(x=Light.Intensity..PAR., y=up, color = Treatment), linewidth=1, show.legend=FALSE)+
geom_ribbon(data=prx, aes(x=Light.Intensity..PAR., ymin=lo, ymax=up, fill=Treatment), alpha = 0.5)
}
p <- p +
scale_color_scico_d(begin=0.9, end=0, palette=a_palette)+
scale_fill_scico_d(begin=0.9, end=0, palette=a_palette)+
labs(title = "Real vs Predicted for MultispeQ logitPS2 ~ Treatment + Light")+
ylab("logit PhiPS2")+
xlab("Light Intensity (PAR)")+
theme_minimal()+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
print(p)
3.1.5 FvP/FmP exploratory
multiPDF_Z(subset(m_data$logitFvPFmP, m_data$logitFvPFmP > 0),
100, c('normal','lognormal','gamma'), a_palette, "FvP/FmP")
multiCDF_Z(subset(m_data$logitFvPFmP, m_data$logitFvPFmP > 0),
100, c('normal','lognormal','gamma'), a_palette, "FvP/FmP")
multiKS_cont(subset(m_data$logitFvPFmP, m_data$logitFvPFmP > 0),
c('normal','lognormal','gamma'))## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Distribution Distance PValue
## 1 Normal 0.068 0.002
## 2 Lognormal 0.057 0.018
## 3 Gamma 0.058 0.015leveneTest(m_data$logitFvPFmP~m_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.595 0.1892
## 718ggplot(data = m_data, aes(x=MinutesFromStart, y = logitFvPFmP,
shape=Treatment, fill=Treatment, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
facet_wrap(~Treatment)+
theme_bw()+
ylab("Logit FvP/FmP")+
xlab("Minutes From Start")+
labs(
title=str_wrap("Logit FvP/FmP By Minutes From Start Across
Inoculation Treatment in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.position = "none",
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.
3.1.6 FvP/FmP analysis
Because FvP/FmP is a unitless ratio from 0:1, we can logit transform it and then treat it as if it were normal.
# Linear model of logit transformed FvPFmP as a function of Treatment with days from germination as a fixed effect
mod_FvP_dfg_s <- summary(mod_FvP_dfg <- (lm(
logitFvPFmP ~ Treatment + DaysFromGermination,
data = m_data)))
print(mod_FvP_dfg_s)##
## Call:
## lm(formula = logitFvPFmP ~ Treatment + DaysFromGermination, data = m_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.271209 -0.065279 -0.004234 0.071809 0.235475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1952892 0.0141777 84.308 < 2e-16 ***
## TreatmentTransplantation 0.0178982 0.0096710 1.851 0.064623 .
## TreatmentGermination -0.0371504 0.0095743 -3.880 0.000114 ***
## TreatmentGerm+Trans -0.0226343 0.0097146 -2.330 0.020087 *
## DaysFromGermination -0.0007802 0.0001833 -4.256 2.36e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09244 on 717 degrees of freedom
## Multiple R-squared: 0.07253, Adjusted R-squared: 0.06736
## F-statistic: 14.02 on 4 and 717 DF, p-value: 5.106e-11## R2 0.06736
# Linear model of logit transformed FvPFmP as a function of Treatment with minutes from start as a fixed effect
mod_FvP_min_s <- summary(mod_FvP_min <- (lm(
logitFvPFmP ~ Treatment + MinutesFromStart,
data = m_data)))
print(mod_FvP_min_s)##
## Call:
## lm(formula = logitFvPFmP ~ Treatment + MinutesFromStart, data = m_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.281106 -0.064749 -0.008379 0.069584 0.231633
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1716386 0.0081069 144.523 < 2e-16 ***
## TreatmentTransplantation 0.0645615 0.0120519 5.357 1.14e-07 ***
## TreatmentGermination -0.0356205 0.0094415 -3.773 0.000175 ***
## TreatmentGerm+Trans 0.0267451 0.0124211 2.153 0.031636 *
## MinutesFromStart -0.0020292 0.0003227 -6.288 5.59e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09112 on 717 degrees of freedom
## Multiple R-squared: 0.0988, Adjusted R-squared: 0.09377
## F-statistic: 19.65 on 4 and 717 DF, p-value: 2.316e-15## R2 0.09377
# Linear model of logit transformed FvPFmP as a function of Treatment with Light Intensity as a fixed effect
mod_FvP_l_s <- summary(mod_FvP_l<- (lm(
logitFvPFmP ~ Treatment + Light.Intensity..PAR.,
data = m_data)))
print(mod_FvP_l_s)##
## Call:
## lm(formula = logitFvPFmP ~ Treatment + Light.Intensity..PAR.,
## data = m_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.250251 -0.070537 -0.007842 0.065663 0.264618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.182e+00 1.127e-02 104.922 < 2e-16 ***
## TreatmentTransplantation 8.900e-03 9.886e-03 0.900 0.368287
## TreatmentGermination -4.478e-02 9.712e-03 -4.611 4.75e-06 ***
## TreatmentGerm+Trans -3.309e-02 9.968e-03 -3.320 0.000947 ***
## Light.Intensity..PAR. -1.398e-04 3.158e-05 -4.428 1.10e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09234 on 717 degrees of freedom
## Multiple R-squared: 0.07442, Adjusted R-squared: 0.06926
## F-statistic: 14.41 on 4 and 717 DF, p-value: 2.521e-11## R2 0.06926
# Linear model of logit transformed FvPFmP as a function of Treatment with Light Intensity AND minutes from start as fixed effects
mod_FvP_lm_s <- summary(mod_FvP_lm<- (lm(
logitFvPFmP ~ Treatment + Light.Intensity..PAR. + MinutesFromStart,
data = m_data)))
print(mod_FvP_lm_s)##
## Call:
## lm(formula = logitFvPFmP ~ Treatment + Light.Intensity..PAR. +
## MinutesFromStart, data = m_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28136 -0.06668 -0.00711 0.06862 0.26022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.218e+00 1.214e-02 100.298 < 2e-16 ***
## TreatmentTransplantation 5.698e-02 1.195e-02 4.770 2.23e-06 ***
## TreatmentGermination -4.383e-02 9.425e-03 -4.650 3.95e-06 ***
## TreatmentGerm+Trans 1.848e-02 1.232e-02 1.500 0.134
## Light.Intensity..PAR. -1.553e-04 3.072e-05 -5.055 5.47e-07 ***
## MinutesFromStart -2.149e-03 3.182e-04 -6.754 2.98e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0896 on 716 degrees of freedom
## Multiple R-squared: 0.1299, Adjusted R-squared: 0.1238
## F-statistic: 21.37 on 5 and 716 DF, p-value: < 2.2e-16## R2 0.069263.2 Fruit
3.2.1 Mass exploratory
# PDF, CDF, and KS test
multiPDF_Z(Fl_data$mass, 100, "all", a_palette, "Mass")
multiCDF_Z(Fl_data$mass, 100, "all", a_palette, "Mass")
multiKS_cont(Fl_data$mass, "all")## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "pexp", rate = x_exp$estimate): ties should not
## be present for the Kolmogorov-Smirnov test## Distribution Distance PValue
## 1 Normal 0.205 0.000
## 2 Lognormal 0.043 0.001
## 3 Gamma 0.104 0.000
## 4 Exponential 0.099 0.000# Tests for homoscedasticity
leveneTest(Fl_data$mass~Fl_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 5.993 0.0004599 ***
## 2080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## plot of mass by treatment with box plots and jittered points
ggplot(data=Fl_data_no_BER, aes(x=Treatment, y=mass, color=Treatment, fill=Treatment))+
geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Tomato Mass (g)")+
xlab("Treatment")+
theme_bw()+
labs(
title="Tomato Mass by Treatment"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot of mass mean per plant by treatment with box plots and jittered points
ggplot(data=Fl_means_byplant, aes(x=Treatment, y=mass_mean, color=Treatment, fill=Treatment))+
geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Mean Tomato Mass (g)")+
xlab("Treatment")+
theme_bw()+
#facet_wrap(~Treatment)+
labs(
title=str_wrap("Mean Tomato Mass Per Plant across Treatments", 50)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.2.2 Mass analysis
According to the PDFs, CDFs, and Kolmogorov-Smirnov tests, mass follows a Gamma or Lognormal distribution. I am unsure of whether to analyze it with each fruit as an observation or with each plant as an observation, where I’m looking at the mean tomato mass per plant instead of the individual tomato masses.
# INDIVIDUAL GAMMA glm
mT_glm_g_s <- summary(mT_glm_g <- glm(mass~Treatment, data=Fl_data, family=Gamma(link="log")))
print(mT_glm_g_s)##
## Call:
## glm(formula = mass ~ Treatment, family = Gamma(link = "log"),
## data = Fl_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.98752 0.04525 88.117 < 2e-16 ***
## TreatmentTransplantation 0.25348 0.06749 3.756 0.000178 ***
## TreatmentGermination 0.06497 0.06692 0.971 0.331704
## TreatmentGerm+Trans -0.04824 0.06490 -0.743 0.457365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.183615)
##
## Null deviance: 1906.7 on 2083 degrees of freedom
## Residual deviance: 1879.7 on 2080 degrees of freedom
## AIC: 21005
##
## Number of Fisher Scoring iterations: 5mT_glm_g_pseudoR2 <- (mT_glm_g_s$null.deviance - mT_glm_g_s$deviance)/mT_glm_g_s$null.deviance
print(mT_glm_g_pseudoR2)## [1] 0.0141688## super high AIC and terrible PR2, so probably not the way to do it.
# GROUPED GAMMA glm
mTg_glm_g_s <- summary(mTg_glm <- glm(mass_mean~Treatment, data=Fl_means_byplant, family=Gamma(link="log")))
print(mTg_glm_g_s)##
## Call:
## glm(formula = mass_mean ~ Treatment, family = Gamma(link = "log"),
## data = Fl_means_byplant)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.81769 0.07606 63.343 <2e-16 ***
## TreatmentTransplantation 0.07786 0.10756 0.724 0.473
## TreatmentGermination 0.07615 0.10756 0.708 0.483
## TreatmentGerm+Trans -0.15292 0.10756 -1.422 0.162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.06941743)
##
## Null deviance: 3.3657 on 47 degrees of freedom
## Residual deviance: 2.9538 on 44 degrees of freedom
## AIC: 472.45
##
## Number of Fisher Scoring iterations: 4mTg_glm_g_pseudoR2 <- (mTg_glm_g_s$null.deviance - mTg_glm_g_s$deviance)/mTg_glm_g_s$null.deviance
print(mTg_glm_g_pseudoR2)## [1] 0.1223773# GROUPED Lognormal glm
mTg_glm_ln_s <- summary(mTg_glm_ln <- glm(mass_mean~Treatment, data=Fl_means_byplant, family=gaussian(link="log")))
print(mTg_glm_ln_s)##
## Call:
## glm(formula = mass_mean ~ Treatment, family = gaussian(link = "log"),
## data = Fl_means_byplant)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.81769 0.07550 63.808 <2e-16 ***
## TreatmentTransplantation 0.07786 0.10286 0.757 0.453
## TreatmentGermination 0.07615 0.10294 0.740 0.463
## TreatmentGerm+Trans -0.15292 0.11593 -1.319 0.194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1046.424)
##
## Null deviance: 52069 on 47 degrees of freedom
## Residual deviance: 46041 on 44 degrees of freedom
## AIC: 475.79
##
## Number of Fisher Scoring iterations: 4mTg_glm_ln_pseudoR2 <- (mTg_glm_ln_s$null.deviance - mTg_glm_ln_s$deviance)/mTg_glm_ln_s$null.deviance
print(mTg_glm_ln_pseudoR2)## [1] 0.1157707# Predictions
## How do I predict these?
# Post Hoc Tests3.2.3 Sugar exploratory
# PDF, CDF, and KS test
multiPDF_Z(Fl_data_no_BER$sugar_avg, 100, "all", a_palette, "Sugar Concentration")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
multiCDF_Z(Fl_data_no_BER$sugar_avg, 100, "all", a_palette, "Sugar Concentration")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
multiKS_cont(Fl_data_no_BER$sugar_avg, "all")## Warning in ks.test.default(x, "pnorm", mean = x_n$estimate[1], sd =
## x_n$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "plnorm", meanlog = x_ln$estimate[1], sdlog =
## x_ln$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced## Warning in ks.test.default(x, "pgamma", shape = x_g$estimate[1], rate =
## x_g$estimate[2]): ties should not be present for the Kolmogorov-Smirnov test## Warning in ks.test.default(x, "pexp", rate = x_exp$estimate): ties should not
## be present for the Kolmogorov-Smirnov test## Distribution Distance PValue
## 1 Normal 0.064 0.012
## 2 Lognormal 0.079 0.001
## 3 Gamma 0.064 0.012
## 4 Exponential 0.436 0.000## plot sugar by treatment as a boxplot with violins and jitter
ggplot(data=Fl_data_no_BER, aes(x=Treatment, y=sugar_avg, color=Treatment, fill=Treatment))+
geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Sugar Concentration (%)")+
xlab("Treatment")+
theme_bw()+
#facet_wrap(~Treatment)+
labs(
title=str_wrap("Tomato Sugar Concentration by Treatment", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot sugar means per plant by treatment as a boxplot with violins and jitter
ggplot(data=Fl_means_byplant, aes(x=Treatment, y=sugar_avg_mean, color=Treatment, fill=Treatment))+
geom_boxplot(width=0.3, alpha=0.5, outliers = FALSE)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Sugar Concentration (%)")+
xlab("Treatment")+
theme_bw()+
#facet_wrap(~Treatment)+
labs(
title=str_wrap("Mean Tomato Sugar Concentration per Plant across Treatments", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## scatter plot of sugar by mass grouped and facet wrapped by treatment
ggplot(data=Fl_data_no_BER, aes(x=mass, y=sugar_avg,
fill=Treatment, shape=Treatment,
color=Treatment))+
geom_point()+
scale_color_manual(values=four_colors)+
scale_fill_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
ylab("Sugar Concentration (%)")+
xlab("Tomato Mass (g)")+
theme_minimal()+
facet_wrap(~Treatment)+
labs(
title="Tomato Sugar Concentration by Mass"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.2.3 Sugar analysis
Sugar appears to have either a Normal or Gamma distribution, so we’ll test both.
# Individual Gamma GLM
mod_sug_g_s <- summary(mod_sug_g <- glm(sugar_avg ~ Treatment, data = Fl_data_no_BER,
family = Gamma(link = "log")))
print(mod_sug_g_s)##
## Call:
## glm(formula = sugar_avg ~ Treatment, family = Gamma(link = "log"),
## data = Fl_data_no_BER)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.86130 0.01719 108.286 < 2e-16 ***
## TreatmentTransplantation 0.03451 0.02393 1.442 0.149722
## TreatmentGermination 0.04273 0.02624 1.628 0.103977
## TreatmentGerm+Trans 0.09020 0.02438 3.699 0.000235 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.04874986)
##
## Null deviance: 34.057 on 627 degrees of freedom
## Residual deviance: 33.375 on 624 degrees of freedom
## AIC: 2311.3
##
## Number of Fisher Scoring iterations: 4# Grouped Gamma glm
mod_sugG_g_s <- summary(mod_sug_g <- glm(sugar_avg_mean ~ Treatment, data = Fl_means_byplant,
family = Gamma(link = "log")))
print(mod_sugG_g_s)##
## Call:
## glm(formula = sugar_avg_mean ~ Treatment, family = Gamma(link = "log"),
## data = Fl_means_byplant)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.85730 0.03448 53.872 <2e-16 ***
## TreatmentTransplantation 0.03351 0.04876 0.687 0.4955
## TreatmentGermination 0.05741 0.04876 1.178 0.2453
## TreatmentGerm+Trans 0.09988 0.04876 2.049 0.0465 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.0142631)
##
## Null deviance: 0.65864 on 47 degrees of freedom
## Residual deviance: 0.59499 on 44 degrees of freedom
## AIC: 117.86
##
## Number of Fisher Scoring iterations: 43.3 Harvest
3.3.1 Harvest exploratory
## scatter plot of ripeness by 'days from harvest to analysis' grouped by treatment (exploratory, to see how delayed analysis may affect sugar content [hopefully it doesn't])
ggplot(data=Fl_data_no_BER, aes(x=DaysFromHarvestToAnalysis, y=ripeness,
fill=Treatment, shape=Treatment,
color=Treatment))+
geom_jitter()+
scale_color_manual(values=four_colors)+
scale_fill_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
ylab("Ripeness")+
xlab("Days from Harvest to Analysis")+
theme_bw()+
#facet_wrap(~Treatment)+
labs(
title="Tomato Ripeness by Days from Harvest to Analysis"
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
## scatter plot of sugar by days from harvest to analysis grouped by treatment (exploratory, to see how delayed analysis may affect sugar content [hopefully it doesn't])
ggplot(data=Fl_data_no_BER, aes(x=DaysFromHarvestToAnalysis, y=sugar_avg,
fill=Treatment, shape=Treatment,
color=Treatment))+
geom_jitter()+
scale_fill_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_color_manual(values=four_colors)+
ylab("Sugar Concentration (%)")+
xlab("Days from Harvest to Analysis")+
theme_bw()+
#facet_wrap(~Treatment)+
labs(
title="Tomato Sugar Concentration by Days from Harvest to Analysis"
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.3.2 Harvest analysis
do tests n shit
## probability of BER for each mass bin
ggplot(data=Fl_data_mb, aes(x=mass_bin, y=pBER*100, fill=mass_bin))+
geom_col()+
ylim(0, 50)+
scale_fill_manual(values=ten_col)+
ylab("BER Probability (%)")+
xlab("Mass Bin (g)")+
theme_bw()+
labs(
title="Tomato BER Probability by Mass Bin"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Plot probability of blossom end rot by treatment as a column graph
ggplot(data=Fl_data_summary, aes(x=Treatment, y=pBER*100, fill=Treatment))+
geom_col()+
scale_fill_manual(values=four_colors)+
ylab("Percent of Fruit w/ BER")+
xlab("Treatment")+
theme_bw()+
labs(
title="Tomato Fruit Blossom End Rot Occurrence by Treatment"
)+
theme(
legend.position="none",
text = element_text(size=30, family="mont", face="bold"),
axis.title = element_text(size=30, family = "mont", face= "bold"),
title = element_text(size=34, family="open", face="bold")
)
## Plot probability of blossom end rot by treatment as a column graph
ggplot(data=Fl_data_summary, aes(x=Treatment, y=fruit_sum, fill=Treatment))+
geom_col()+
scale_fill_manual(values=four_colors)+
ylab("Fruit Count")+
xlab("Treatment")+
theme_bw()+
labs(
title="Total Fruit Count by Inoculation Treatment"
)+
theme(
legend.position="none",
text = element_text(size=30, family="mont", face="bold"),
axis.title = element_text(size=30, family = "mont", face= "bold"),
title = element_text(size=34, family="open", face="bold")
)
## Plot blossom end rot by treatment as a stacked bar chart
ggplot(data=Fl_data, aes(x=Treatment, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
geom_bar(position="stack", stat="identity")+
scale_fill_manual(values=two_colors)+
ylab("Fruit Count")+
xlab("Treatment")+
theme_bw()+
guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
labs(
title="Tomato Fruit and BER Count by Treatment"
)+
theme(
legend.position="right",
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## Plot blossom end rot by mass bin as a stacked bar chart
ggplot(data=Fl_data, aes(x=mass_bin, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
geom_bar(position="stack", stat="identity")+
scale_fill_manual(labels=c("No", "Yes"), values=two_colors)+
guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
ylab("Fruit Count")+
xlab("Mass Bin")+
theme_bw()+
labs(
title="Tomato Fruit and BER Count by Mass Bin"
)+
theme(
legend.position="inside",
legend.position.inside=c(.8,.7),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.text.x = element_text(size=20, family = "mont", angle=45,
lineheight=1, hjust=1, vjust=1.3),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## Plot blossom end rot by plant as a stacked bar chart
ggplot(data=Fl_data, aes(x=plant, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
geom_bar(position="stack", stat="identity")+
scale_fill_manual(labels=c("No", "Yes"), values=two_colors)+
guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
ylab("Fruit Count")+
xlab("Plant")+
facet_wrap(~Treatment)+
theme_bw()+
labs(
title="Tomato Fruit and BER Count by Plant"
)+
theme(
legend.position="right",
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.text = element_text(size=16, family = "mont",
lineheight=1),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## fruit mass sum by treatment
ggplot(data=Fl_data_summary, aes(x=Treatment, y=mass_sum/1000,
fill=Treatment))+
geom_col()+
scale_fill_manual(values=four_colors)+
ylab("Total Fruit Mass (kg)")+
xlab("Treatment")+
theme_bw()+
labs(
title="Total Tomato Fruit Harvest across Inoculation Treatments"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## fruit mass sum by plant
ggplot(data=Fl_sum_byplant, aes(x=plant, y=mass_sum/1000,
fill=Treatment))+
geom_col()+
facet_wrap(~Treatment)+
scale_fill_manual(values=four_colors)+
scale_x_discrete(name="Plant", breaks=c("3","6","9"))+
ylab("Total Fruit Mass (kg)")+
xlab("Plant")+
theme_bw()+
labs(
title="Total Tomato Fruit Harvest across Individual Plants"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot of fruit mass by treatment
ggplot(data=Fl_sum_byplant, aes(x=Treatment, y=mass_sum/1000, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5, outliers = FALSE)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Mass Sum (kg)")+
xlab("Treatment")+
theme_bw()+
labs(
title="Total Fruit Mass per Plant across Inoculation Treatments"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot of fruit count by treatment
ggplot(data=Fl_sum_byplant, aes(x=Treatment, y=fruit_sum, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5, outliers = FALSE)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Fruit Count")+
xlab("Treatment")+
theme_bw()+
labs(
title="Tomato Plant Harvested Fruit Count Across Inoculation Treatments"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
Fl_avg_sums <- Fl_sum_byplant %>%
group_by(Treatment) %>%
summarize_at(vars(mass_sum), list(mean=mean, sd=sd))%>%
mutate(mean=mean/1000,
sd=sd/1000)
## mass per treatment columns with error bars
ggplot(data=Fl_avg_sums, aes(x=Treatment, y=mean, fill=Treatment))+
geom_col()+
geom_errorbar( aes(ymin=mean-sd, ymax=mean+sd), width=0.4)+
scale_fill_manual(values=four_colors)+
ylab("Total Mass (kg)")+
xlab("Treatment")+
labs(
title="Total Tomato Fruit Mass (kg) Across Inoculation Treatments"
)+
theme_bw()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)